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Abstract 

In the present review we survey the properties of a transcendental function of the Wright 
type, nowadays known as M-Wright function, entering as a probability density in a 
relevant class of self-similar stochastic processes that we generally refer to as time- 
fractional diffusion processes. Indeed, the master equations governing these processes 
generalize the standard diffusion equation by means of time-integral operators interpreted 
as derivatives of fractional order. When these generalized diffusion processes are properly 
characterized with stationary increments, the M- Wright function is shown to play the 
same key role as the Gaussian density in the standard and fractional Brownian motions. 
Furthermore, these processes provide stochastic models suitable for describing phenomena 
of anomalous diffusion of both slow and fast type. 

1 Introduction 

By time-fractional diffusion processes we mean certain diffusion- like phenomena governed 
by master equations containing fractional derivatives in time whose fundamental solution 
can be interpreted as a probability density function {pdf) in space evolving in time. It 
is well known that for the most elementary diffusion process, the Brownian motion, the 
master equation is the standard linear diffusion equation whose fundamental solution is the 
Gaussian density with a spatial variance growing linearly in time. In such case we speak 
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about normal diffusion, reserving the term anomalous diffusion when the variance grows 
differently. A number of stochastic models for explaining anomalous diffusion have been 
introduced in literature, among them wc like to quote the fractional Brownian motion, 
see e.g. [50, 70], the Continuous Time Random Walk, see e.g. [25, 51, 53, 63], the Levy 
flights, see e.g. [11], the Schneider grey Brownian motion, see [64, 65], and, more generally, 
random walk models based on evolution equations of single and distributed fractional 
order in time and/or space, see e.g. [7, 8, 9], [23, 24], [33, 34], [77, 78]. 

In this survey paper we focus our attention on modifications of the standard diffusion 
equation, where the time can be stretched by a power law (t — )■ < a < 2) and the first- 
order time derivative can be replaced by a derivative of non- integer order /3 (0 < /3 < 1). 
In these cases of generalized diffusion processes the corresponding fundamental solution 
still keeps the meaning of a spatial pdf evolving in time and is expressed in terms of 
a special function of the Wright type that reduces to the Gaussian when P — 1. This 
transcendental function, nowadays known as M- Wright function, will be shown to play 
a fundamental role for a general class of self-similar stochastic processes with stationary 
increments, which provide stochastic models for anomalous diffusion, as recently shown 
by Mura et al. [55, 56, 57, 58]. 

In Section 2 we provide the reader with the essential notions and notations concerning the 
integral transforms and fractional calculus, which are necessary in the rest of the paper. 
In Section 3 we introduce in the complex plane C the series and integral representations of 
the general Wright function denoted by Wx^i^{z) and of the two related auxiliary functions 
F,,{z), M,y{z), which depend on a single parameter. In Section 4 we consider our auxiliary 
functions in real domain pointing out their main properties involving their integrals and 
their asymptotic representations. Mostly, we restrict our attention to the second auxiliary 
function, that we call M- Wright function, when its variable is in IR + or in all of IR but 
extended in symmetric way. We derive a fundamental formula for the absolute moments 
of this function in IR which allows us to obtain its Laplace and Fourier transforms. 
In Section 5 we consider some types of generalized diffusion equations containing time 
partial derivatives of fractional order and we express their fundamental solutions in 
terms of the M- Wright functions evolving in time with a given self-similarity law. In 
Section 6 we stress how the M- Wright function emerges as a natural generalization of the 
Gaussian probability density for a class of self-similar stochastic processes with stationary 
increments, depending on two parameters {a, f3). These processes are defined in a unique 
way by requiring the determination of any multi-point probability distribution and include 
the well-known standard and fractional Brownian motion. We refer to this class as the 
generalized grey Brownian motion (ggBm), because it generalizes the grey Brownian 
motion {gBm) introduced by Schneider [64, 65]. Finally, a short concluding discussion is 
drawn. In Appendix A we derive the fundamental solution of the time-fractional diffusion 
equation. In Appendix B we outline the relevance of the M- Wright function in time- 
fractional drift processes entering as subordinators in time-fractional diffusion. 
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2 Notions and Notations 

Integral transforms pairs. 

In our analysis we will make extensive use of integral transforms of Laplace, Fourier and 
Mellin type so we first introduce our notation for the corresponding transform pairs. We 
do not point out the conditions of validity and the main rules, since they are given in any 
textbook on advanced mathematics. 

Let 

f{s)=/:{firy,r^s}= e-^V(r)dr, (2.1) 

Jo 

be the Laplace transform of a sufficiently well-behaved function /(r) with r e IR"^, s e C, 
and let 

/(r) = {/(.); s ^ r } = 2^ ^ e f{s) ds , (2.2) 

be the inverse Laplace transform, where Br denotes the so-called Bromwich path, a 
straight line parallel to the imaginary axis in the complex s-plane. Denoting by •<-> the 
justaposition of the original function /(r) with its Laplace transform f{s), the Laplace 
transform pair reads 

fir) A fis). (2.3) 

Let 

/ + 0O 
e+"^^f{x)dx, (2.4) 
-oo 

be the Fourier transform of a sufficiently well-behaved function f{x) with a; e IR, k e IR, 
and let 

f{x)^:F-'[f{Ky,K^x]^^ e-'^^f{K)dK, (2.5) 

be the inverse Fourier transform. Denoting by •H- the justaposition of the original function 
f{x) with its Fourier transform /(«), the Fourier transform pair reads 

fix) 4 fin). (2.6) 

Let 

/•oo 

f*{s)^M{fir)-r^s}^ r'-^fir)dr, (2.7) 

Jo 

be the Mellin transform of a sufficiently well-behaved function /(r) with r e IR"*", s e (C, 
and let 

fir) = {f*is);s ^r}^^J^r-' f*is) ds , (2.8) 
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be the inverse MeUin transform. Denoting by the justaposition of the original function 
/(r) with its MeUin transform the MeUin transform pair reads 

m ^ ris). (2.9) 



Essentials of fractional calculus with support in IR"*" . 

Fractional calculus is the branch of mathematical analysis that deals with pseudo- 
differential operators that extend the standard notions of integrals and derivatives to any 
positive non-integer order. The term fractional is kept only for historical reasons. Let 
us restrict our attention to sufficiently well-behaved functions f{t) with support in IR"*". 
Two main approaches exist in the literature of fractional calculus to define the operator 
of derivative of non integer order for these functions, referred to Riemann-Liouville and 
to Caputo. Both approaches are related to the so-called Riemann-Liouville fractional 
integral defined for any order > as 

Jtfit) ^ - rr-'f{r) dr . (2.10) 
We note the convention = I (Identity) and the semigroup property 

j'^ ^ ji; ^ /.>o,i/>o. (2.11) 

The fractional derivative of order > in the Riemann-Liouville sense is defined as the 
operator which is the left inverse of the Riemann-Liouville integral of order ji (in 
analogy with the ordinary derivative), that is 

D^J.^^I, //>0. (2.12) 

If m denotes the positive integer such that m— l</x<m,we recognize from Eqs. (2.11) 
and (2.12): f{t) := JT~^ f{t) , hence 



dt"" 
d"" 



1 /■* f{T)dT 



dt"^ 

For completeness we define D? — I. 



T{m - n) Jo {t - r)^+i-"^ 
f{t) , // = m 



m — 1 < II < m, 

(2.13) 



On the other hand, the fractional derivative of order > in the Caputo sense is defined 
as the operator ^D^ such that f{t) := J^'" D"^ f{t) , hence 



/* 

^0 



/("*)(t) dr 

rn — u) In I 

d" 



m — 1 < /I < m, 



^D'l fit) = { - /^) Jo (t - t)^+i- ' ^ ' (2.14) 

;f{t) , fi^m. 



dt'' 
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We note the different behavior of the two derivatives in the hmit — > (m — 1)+. In fact, 



— >■ (m — 1)" 



fit) 



D^f{t) ^ Dr Ji fit) = a 

.D^fit) ^ DT fit) = Df-'^ fit) - Dr-"f iO+) , 



,(m-l) 



(2.15) 



where the hmit for i — > 0"^ is taken after the operation of derivation. 

Furthermore, recaUing the Riemann-Liouville fractional integral and derivative of the 
power law for i > 0, 



r(7 + i) 

r(7 + i + /^) 
r(7 + i) 



A(>0, 7> -1, 



(2.16) 



r(7 + 1 - fi) 

we find the relationship between the two types of fractional derivative. 



f(t)-j:^j''\^ 



k=0 



^DUit). 



(2.17) 



We note that the Caputo definition for the fractional derivative incorporates the initial 
values of the function and of its integer derivatives of lower order. The subtraction of 
the Taylor polynomial of degree m — 1 at f = 0"^ from fit) is a sort of regularization of 
the fractional derivative. In particular, according to this definition, the relevant property 
that the derivative of a constant is zero is preserved for the fractional derivative. 

Let us finally point out the rules for the Laplace transform with respect to the fractional 
integral and the two fractional derivatives. These rules are expected to properly generalize 
the well-known rules for standard integrals and derivatives. 

For the Riemann-Liouville fractional integral we have 



c{jn{t)-,t^s} 



fXf) 



For the Caputo fractional derivative we consequently get 



m—1 



£ {.D^ fit);t ^s}^s'^fis)-J2 /^'nO+) ,m -!<!,< 



m , 



(2.18) 



(2.19) 



fe=0 



where /^'^^(0+) := lim f^''\t). The corresponding rule for the Riemann-Liouville 
fractional derivative is more cumbersome and it reads 



m—1 



C{D^fit);t^s}^s> 



(m-fj.) 



k=0 



f{0' 



^m—l—k 



, m-1 < II <m, (2.20) 
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where the hmit for t — > 0+ is understood to be taken after the operations of fractional 
integration and derivation. As soon as all the limiting values /'^^•*(0+) are finite and 
m — 1 < II < m, formula (2.20) for the Riemann-Liouville derivative simplifies into 

C{D'^ f{t);t^ s} = s^fis), m-l< ii<m. (2.21) 

In the special case f^^\Q'^) = for A; = 0, 1, m — 1, we recover the identity between the two 
fractional derivatives. The Laplace transform rule (2.19) was practically the key result of 
Caputo [5, 6] in defining his generalized derivative in the late sixties. The two fractional 
derivatives have been well discussed in the 1997 survey paper by Gorenflo and Mainardi 
[21], see also [42], and in the 1999 book by Podlubny [59]. In these references the Authors 
have pointed out their preference for the Caputo derivative in physical applications where 
initial conditions are usually expressed in terms of finite derivatives of integer order. 

For further reading on the theory and applications of fractional calculus we recommend 
the recent treatise by Kilbas et al. [29] . 



3 The functions of the Wright type 



The general Wright function. 

The Wright function, that we denote by Wx,^j.X^)i named in honour of E. Maitland 
Wright, the eminent British mathematician, who introduced and investigated this function 

in a series of notes starting from 1933 in the framework of the asymptotic theory of 
partitions, see [73, 74, 75]. The function is defined by the series representation, convergent 
in the whole 2;- complex plane, 

00 n 

WxM■■=Y.^¥^r—-^^ A>-l,//eC. (3.1) 
^ n\ T[\n + //) 

Originally, Wright assumed A > 0, and, only in 1940 [76], he considered — 1 < A < 0. We 
note that in Chapter 18 of Vol. 3 of the handbook of the Bateman Project [12], devoted 
to Miscellaneous Functions, presumably for a misprint, the parameter A of the Wright 
function is restricted to be non negative. When necessary, we propose to distinguish the 
Wright functions in two kinds according to A > {first kind) and — 1 < A < {second 
kind). 

For more details on Wright functions the reader can consult e.g. [19, 20, 28, 30, 39, 69, 
71, 72] and references therein. 

The integral representation of the Wright function reads 
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where Ha denotes the Hankel path. We remind that the Hankel path is a loop that starts 
from — oo along the lower side of the negative real axis, encircles the circular area around 
the origin with radius e — ?> in the positive sense, and ends at — oo along the upper side of 
the negative real axis. The equivalence of the series and integral representations is easily 
proved using Hankel formula for the Gamma function 



1 



r(C) 



J Ha 



and performing a term-by-term integration. In fact. 



2m 



Ha 



-A da 



1 

27ri 



Ha 



E 

,n=0 



a 



-An 



da 



E 

n=0 



Z 



1 
27ri 



Ha 



E 

n=0 



n\ T[\n -\- /x] 



It is possible to prove that the Wright function is entire of order 1/(1 + A) , hence it is of 
exponential type only if A > (which corresponds to Wright functions of the first kind). 
The case A = is trivial since Wq^^^z) = e^/r(/i) , provided that 7^ 0, —1, —2, .... 



The auxiliary functions of the Wright type. 

Mainardi, in his first analysis of the time- fractional diffusion equation [36, 48], aware of 
the Bateman handbook [12], but not yet of the 1940 paper by Wright [76], introduced the 
two (Wright-type) entire auxiliary functions, 

F,{z):=W_,,o{-z), 0<z/<l, (3.3) 

and 

M,{z) := W.,,i.,{-z) , < i/ < 1 , (3.4) 

inter-related through 

F,{z) = ly z M,{z) . (3.5) 

As a matter of fact, functions Fi,{z) and M^(z) are particular cases of the Wright function 
of the second kind Wx^i^{z) by setting A = —u and /i — or /i — 1, respectively. 

Hereafter, we provide the series and integral representations of the two auxiliary functions 
derived from the general formulas (3.1) and (3.2), respectively. 

The series representations for the auxiliary functions read 

FJz) := y , = - y Vivn + 1) sini-Kvu) , (3.6) 

n=l ^ ' n=l 
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and 

oo / \^ -J oo / \n—l 

MJz) := Y ^ = -Y 7^^^^ r(z/n) sin(7rz/n) . (3.7) 

^ n\T[-vn + il - v)] t: ^ in - 1)1^ ' ^ ' ^ ' 

n=0 ^ ^ n=l ^ ' 

The second series representations in Eqs. (3.6)-(3.7) have been obtained by using the 
reflection formula for the Gamma function T{C,) r(l — Q — n/ sin ttC- 

As an exercise, the reader can directly prove that the radius of convergence of the power 
series in (3.6)-(3.7) is infinite for < z/ < 1 without being aware of Wright's results, as it 
was shown independently by Mainardi [36], see also [59]. 

Furthermore, we have F,y{0) = and M^(0) = l/r(l — u). We note that relation (3.5) 
between the two auxiliary functions can be easily deduced from (3.6)-(3.7), by using the 
basic property of the Gamma function r(^ + 1) = (^r(^). 

The integral representations for the auxiliary functions read 

F,{z):=^ f e'^-^'^'^da, (3.8) 



27ri 

M, 



Ha 



We note that relation (3.5) can be obtained also from (3.8)-(3.9) with an integration by 
parts. In fact, 

M^{z) ^ f e^-'^' ^ = [ (-^^e-'^'] da 

J Ha (J^-'' J Ha \ da 



VZ J Ha ^Z 

The equivalence of the scries and integral representations is easily proved by using the 
Hankel formula for the Gamma function and performing a term-by-term integration. 



Special cases. 

Explicit expressions of F^i^z) and M,^{z) in terms of known functions are expected for 
some particular values of u. Mainardi and Tomirotti [48] have shown that for u — 1/q, 
where q > 2 is a positive integer, the auxiliary functions can be expressed as a sum of 
simpler {q — 1) entire functions. In the particular cases q = 2 and g = 3 we find 

1 °° / 1 \ r^"^ 1 
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and 

M,/S{Z) -Ym^XyrnJ^-' Wmj^XVrnWn+iy. (3.11) 

= 3'/'Ai (^/3V3), 
where Ai denotes the Airy function. 

Furthermore, it can be proved that Mi/g{z) satisfies the differential equation of order q — 1 

(-1)1 

^ Mi/,iz) + ^--^zMy.iz) = , (3.12) 
subjected to the q — 1 initial conditions at 2; = 0, derived from (3.7), 



O) = i (1^1 ^y^^ = m + 1)/?] sin[7r {h + l)/q] , (3.13) 



with h = 0, 1, ... g — 2. We note that, for g > 4, Eq. (3.12) is akin to the hyper- 
Airy differential equation of order q — 1, see e.g. [3]. Consequently, the auxiliary 
function Mi,{z) could be considered as a sort of generalized hyper-Airy function. However, 
in view of further applications in stochastic processes, we prefer to consider it as a 
natural (fractional) generalization of the Gaussian function, similarly as the Mittag-Lefiler 
function is known to be the natural (fractional) generalization of the exponential function. 
To stress the relevance of the auxihary function My{z), it was also suggested the special 
name M- Wright function, a terminology that has been followed in literature to some 
extent^. 

Moreover, the analysis of the limiting cases u — and u — 1 requires special attention. 
For 1/ — we easily recognize from the series representations (3.6)-(3.7): 

Foiz)=0, Mo{z)=e-\ 

The limiting case u = 1 is singular for both the auxiliary functions as expected from the 
definition of the general Wright function when A = -i^ = —1. Later we will deal with 
this singular case for the M- Wright function when the variable is real and positive. 



^Some authors including Podlubny [59], Gorenflo et al. [19, 20], Hanyga [26], Balescu [2], Chechkin et 
al. [9], Germane et al. [16], Kiryakova [31, 32] refer to the Af- Wright function as the Mainardi function. 
It was Professor Stankovic, during the presentation of the paper by Mainardi and Tomirotti [48] at the 
Conference Transform Methods and Special Functions, Sofia 1994, who informed Mainardi, being aware 
only of the Bateman Handbook [12], that the extension for —1 < A < had been already made just by 
Wright himself in 1940 [76], following his previous papers published in the thirties. Mainardi, in the paper 
[43] devoted to the 80-th birthday of Prof. Stankovic, used the occasion to renew his personal gratitude 
to Prof. Stankovic for this earlier information that led him to study the original papers by Wright and 
work (also in collaboration) on the functions of the Wright type for further applications, see e.g. [19, 20] 
and [451. 
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4 Properties and plots of the auxiliary Wright 
functions in real domain 

Let us state some relevant properties of the auxiliary Wright functions, with special 
attention to the M^, function in view of its role in time-fractional diffusion processes. 

Exponential Laplace transforms. 

We start with the Laplace transform pairs involving exponentials in the Laplace domain. 
These were derived by Mainardi in his earlier analysis of the time fractional diffusion 
equation, see e.g. [36], [37], 

- {l/r^ -^M, A e-^\ Q<u<l, (4.1) 

-F,(1/0 = 1m,(1/0 4 0<u<l. (4.2) 

u r s 

We note that the inversion of the Laplace transform of the exponential exp {—s") is 

relevant since it yields for any u G (0, 1) the unilateral extremal stable densities in 

probability theory, denoted by L~'^{r) in [44]. As a consequence, the non-negativity 

of both the auxiliary Wright functions when their argument is positive is proved by the 

Bernstein theorem^. The Laplace transform pair in (4.1) has a long history starting from 

a formal result by Humbert [27] in 1945, of which Pollard [61] provided a rigorous proof 

one year later. Then, in 1959 Mikusihski [54] derived a similar result on the basis of his 

theory of operational calculus. In 1975, albeit unaware of the previous results, Buchen 

and Mainardi [4] derived the result in a formal way. We note that all the above authors 

were not informed about the Wright functions. To our actual knowledge the former author 

who derived the Laplace transforms pairs (4.1)-(4.2) in terms of Wright functions of the 

second kind was Stankovic in 1970, see [69]. 

Hereafter we would like to provide two independent proofs of (4.1) carrying out the 
inversion of exp (— s^) , either by the complex Bromwich integral formula following [36], or 
by the formal series method following [4]. Similarly we can act for the Laplace transform 
pair (4.2). For the complex integral approach we deform the Bromwich path Br into the 
Hankel path Ha, that is equivalent to the original path, and we set a — sr. Recalling the 
integral representation (3.8) for the Fi, function and Eq. (3.5), we get 

>C-Mexp(-s^);s^r] = -^ [ e''' ' ''^ ds ^ f e"" ' (""M'^ da 

= ^F.(l/0 = -^M.(l/0. 
^We refer to Feller's treatise [13] for Laplace transforms, stable densities and Bernstein theorem. 
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Expanding in power series the Laplace transform and inverting term by term, we formally 
get 

£-1 [exp i-sn] = y ^—^ [s^ = y — - 

^ ^ ^ ^ n\ ^ ^ ^ n\ Ti-vn) 

n=0 n=l ^ ' 

= ^F.(l/0 = -^M.(l/0, 

where now we have used the series representation (3.6) for the function F^, along with the 
relationship formula (3.5). 



Asymptotic representation for large argument. 

Let us point out the asymptotic behaviour of the function Mi,(r) when r — )■ cxd. Choosing 
as a variable r/v rather than r, the computation of the desired asymptotic representation 
by the saddle-point approximation is straightforward. Mainardi and Tomirotti [48] have 
obtained 

M,{r/v) ~ a{v) " V2)/(l - i^) gxp \-h{v) r^l ~ , 

(4.3) 

aiv) = , > , biv) = > . 



The above evaluation is consistent with the first term in the asymptotic series expansion 
provided by Wright with a cumbersome and formal procedure for his general function VTa,/* 
when — 1 < A < 0, see [76]. In 1999 Wong and Zhao have derived asymptotic expansions 
of the Wright functions of the first and second kind in the whole complex plane following 
a new method for smoothing Stokes' discontinuities, see [71, 72], respectively. 

We note that, for v — 1/2 Eq. (4.3) provides the exact result consistent with (3.10), 

Mv2(2r) = ^e-^' ^ My^{r) = -^e'"''!^ . (4.4) 

We also note that in the limit p ^ 1~ the function M^[r) tends to the Dirac generalized 
function 5 (r — 1), as can be recognized also from the Laplace transform pair (4.1). 



Absolute moments. 

Prom the above considerations we recognize that, for the M-Wright functions, the 
following rule for absolute moments in IR+ holds 

rr'M,{r)dr^^j^-^, S > -1 , < u < 1 . (4.5) 

Jo r(i/(^ + 1) 
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In order to derive this fundamental result, we proceed as follows on the basis of the integral 
representation (3.9): 



rMJr)dr = 





1 

2wi 



1 
27ri 



e 

Ha 



da 
^1—1/ 



c 

Ha 



dr 
da 



a 



l-u 



r(^+i) 

27ri 



Ha 



^us+i Y{v5 + 1) ■ 



da 



Above we have legitimized the exchange between integrals and used the identity 



along with the Hankel formula of the Gamma function. Analogously, we can compute all 
the moments of Fi,(r) in IR"*". 



The Laplace transform of the M- Wright function. 

Let the Mittag-Leffler function be defined in the complex plane for any > by the 
following series and integral representation, see e.g. [12, 41], 

^r(i/n+l) 2m J Ha -z 

Such function is entire of order 1/a for a > and reduces to the function exp {z) for 
z/ > and to 1/(1 — 2;) for z/ = 0. We recall that the Mittag-Lcffler function for z/ > 
plays fundamental roles in applications of fractional calculus like fractional relaxation and 
fractional oscillation, see e.g. [1], [21], [42], [40], so that it could be referred as the Queen 
function of fractional calculus^. 

We now point out that the M- Wright function is related to the Mittag-LefHer function 
through the following Laplace transform pair, 

M^(r) h E^{-s) , < z/ < 1 . (4.7) 

For the reader's convenience we provide a simple proof of (4.7) by using two different 
approaches. AW assume that the exchanges between integrals and series are legitimate in 

^Recently, numerical routines for functions of Mittag-Leffler type have been provided e.g. by Freed 
et al. [14], Gorenflo et al. [18] (with MATHEMATICA), Podlubny [60] (with MATLAB), Seybold and 
Hilfer [67]. 
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view of the analyticity properties of the involved functions. In the first approach we use 
the integral representations of the two functions obtaining 



^0 



, —sr 



MJr)dr 



,(T — rcr 



da 



.l-v 



1 1'°° 
27rz Jo 

= — I e^a^-^ [re-^(« + Odr 

27rZ J Ha Uq 



da 



(4.8) 



-—I 



a „v — 1 



da — Ei,{—s) . 



In the second approach we develop in series the exponential kernel of the Laplace transform 
and we use the expression (4.5) for the absolute moments of the M- Wright function 
arriving to the following series representation of the Mittag-Leffler function, 



poo °° / \n roo 

/ e-^''M,{r)dr = / r''M,{r)dr 

Jo „_n ^' Jo 



n=0 

oo 

E 

n=0 



-s)" T(n + 1) 



E 

n=0 



(4.9) 



Tivn+l] 



We note that the transformation term by term of the series expansion of the M- Wright 
function is not legitimate because the function is not of exponential order, see [10]. 
However, this procedure yields the formal asymptotic expansion of the Mittag-Leffler 
function Ei,{—s) as s — >■ oo in a sector around the positive real axis, see e.g. [12, 41], that 
is 

-^^(-r)"rfr ^ (-ir 1 



oo roo 

E Jo ^ 

n=0 ^ ^ 



E 

n=0 

oo 



r(-z/n + 1 - /y) 



r(-i/m + 1) s*^ 



~ Ei,[—s) , s — >■ oo . 



The Fourier transform of the symmetric M- Wright function. 

The M- Wright function, extended on the negative real axis as an even function, is related 
to the Mittag-Leffler function through the following Fourier transform pair 



M,{\x\) A 2E2,{-k'), 0<i/< 1. 
Below, we prove the equivalent formula 

poo 

/ cos(Kr) M„{r) dr = E2^{-k'^) . 
Jo 



(4.10) 



(4.11) 
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For the prove it is sufficient to develop in series the cosine function and use formula (4.5) 
for the absolute moments of the M- Wright function, 

/•oo 2n /-oo 

/ cosUr) MJr) dr =V(-1)»^ — - / r'"MJr)dr 
Jo ti (^^yJ" (4.12) 



oo 



n=0 ^ ^ 



The Mellin transform of the M- Wright function. 

It is straightforward to derive the Mellin transform of the M- Wright function using result 
(4.5) for the absolute moments of the M- Wright function. In fact, setting S — s — 1 in 
(4.5), by analytic continuation it follows 



Plots of the symmetric M- Wright function. 

It is instructive to show the plots of the (symmetric) M- Wright function on the real axis 
for some rational values of the parameter u. In order to have more insight of the effect 
of the parameter itself on the behaviour close to and far from the origin, we adopt both 
linear and logarithmic scale for the ordinates. 

In Figs. 1 and 2 we compare the plots of the M;^(x)- Wright functions in —5 < a; < 5 for 
some rational values of u in the ranges u G [0, 1/2] and u G [1/2, 1], respectively. In Fig. 
1 we see the transition from exp (— |a;|) for z/ = to 1/y^exp (— x^) for u = 1/2, whereas 
in Fig. 2 we see the transition from l/^/^Texp (— x^) to the delta functions 5{x ±1) for 
u — 1. Because of the two symmetrical humps for 1/2 < i/ < 1, the function appears 
bi-modal with the characteristic shape of the capital letter M. 

In plotting M,y{x) at fixed u for sufficiently large x the asymptotic representation (4.3)- 
(4.4) is useful since, as x increases, the numerical convergence of the series in (3.7) 
decreases up to being completely inefficient: henceforth, the matching between the series 
and the asymptotic representation is relevant and followed by Mainardi and associates, 
see e.g. [38, 39, 44, 45]. However, as — >■ 1~, the plotting remains a very difficult task 
because of the high peak arising around x — ±1. For this we refer the reader to the 1997 
paper by Mainardi and Tomirotti [49], where a variant of the saddle point method has 
been successfully used to properly depict the transition to the delta functions 6{x ± 1) 
as i> approaches 1. For the numerical point of view we like to highlight the recent paper 
by Luchko [35], where algorithms are provided for computation of the Wright function on 
the real axis with prescribed accuracy. 
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Figure 1: Plots of the symmetric Mj,- Wright function with v = 0,1/8,1/4,3/8,1/2 for 
— 5 < X < 5; left: linear scale, right: logarithmic scale. 




X X 



Figure 2: Plots of the symmetric M- Wright function with u = 1/2, 5/8, 3/4, 1 for 
—5 < X < 5; left: linear scale; right: logarithmic scale. 

The M- Wright function in two variables. 

In view of the time-fractional diffusion processes that will be considered in the next 
Sections, it is worthwhile to introduce the function in two variables 

M^{x,t) ■=r'' M^ixr") , o<z/<i, x,teWi^, (4.14) 

which defines a spatial probability density in x evolving in time t with self-similarity 
exponent H = v. Of course for x G IR we have to consider the symmetric version 
obtained from (4.14) multiplying by 1/2 and replacing x by 

Hereafter we provide a list of the main properties of this function, which can be derived 
from Laplace and Fourier transforms of the corresponding M- Wright function in one 
variable. 

From Eq. (4.2) we derive the Laplace transform of lMi,(x,t) with respect to t G IR^, 

£{lM^(a;,t);t ^ s} = s^-^e~^^^ (4.15) 
From Eq. (4.6) we derive the Laplace transform of lMj^(x,t) with respect to x G 1R + , 



C {lM^(a;, t)]x ^ s} = {.st") 



(4.16) 
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Prom Eq. (4.10) we derive the Fourier transform of IM,^(|a;|, t) with respect to x e IR, 

{M^{\x\,t);x k} = 2E2^ {-kH") . (4.17) 

Moreover, using the MeUin transform, Mainardi et al. [46] derived the following integral 
formula, 

poo 

M^{x,t)= Mx{x,T)M^{T,t)dT, u = \ii. (4.18) 

Special cases of the M- Wright function are simply derived for z/ = 1/2 and u = 1/3 
from the corresponding ones in the complex domain, see Eqs. (3.10)-(3.11). We devote 
particular attention to the case i/ = 1/2 for which we get from (4.4) the Gaussian density 
in IR, 

lM.,.(|.|,() = ^e-V(«). (4.19) 
For the limiting case = 1 we obtain 

^M,{\x\,t) = ^ [Six -t) + Six + t)] . (4.20) 



5 Fractional diffusion equations 

Let us now consider a variety of diffusion-like equations starting from the standard 
diffusion equation whose fundamental solutions are expressed in terms of the M- Wright 
function depending on space and time variables. The two variables, however, turn out to 
be related through a self-similarity condition. 



The standard diffusion equation. 

The standard diffusion equation for the field uix, t) with initial condition uix, 0) = uoix) 
is 

— = -oo<x<oo,t>0, (5.1) 

where Ki is a suitable diffusion coefficient of dimensions [Ki\ = [L]^[T]~^ = crn?/sec. 
This initial-boundary value problem can be easily shown to be equivalent to the Volterra 
integral equation 

uix, t) = Uoix) + ' dr . (5.2) 

It is well known that the fundamental solution (usually refereed as the Green function), 
which is the solution corresponding to uoix) = Six), is the Gaussian probability density 
evolving in time with variance (mean square displacement) proportional to time. In our 
notation we hve: 
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/ + 00 
Qi{x,t) dx = 2Kit. 
-oo 



(5.4) 



This variance law characterizes the process of normal diffusion as it emerges from 
Einstein's approach to Brownian motion {Bm), see e.g. [68]. 

In view of future developments, we rewrite the Green function in terms of the M- Wright 
function by recalling Eq. (3.10), that is. 



Prom the self-similarity of the Green function in (5.3) or (5.5) we are led to write 



Gi{x,t) 



/Kit" 



\x\ 



/Kit" 



,1 



(5.5) 



(5.6) 



where H = 1/2 is the similarity (or Hurst) exponent and ^ = \x\/ [\/Klt^^'^) acts as 
the similarity variable. We refer to the one- variable function Qi as the reduced Green 
function. 



The stretched-time standeird diffusion equation. 

Let us now stretch the time variable in Eq. (5.1) by replacing t with where < a < 2. 
We have 

-oo<x<+oo, i>0, (5.7) 

where K^ is a sort of stretched diffusion coefficient of dimensions [Ka\ = [L]^[T]~" = 
cm^/sec"'. It is easy to recognize that such equation is akin to the standard diffusion 
equation but with a diffusion coefficient depending on time, Ki{t) — at°'~^ K^- In fact, 
using the rule 

d 1 d 



dt<^ at<^-^ dt ' 

we have 

j.a—1 



du „, 1 ^, d'^u 



— = af^-^ KaTT^, -oo<x < +00 , t>0. (5.8) 
at ox^ 

The integral form corresponding to Eqs. (5.7)-(5.8) reads 

u{x, t) = uo{x) + aK^ ^!^llr"-^ dr . (5.9) 



The corresponding fundamental solution is the stretched-time Gaussian 

Ga{x,t) = J_ .-^V(4i^X) = i_i_M,/, ( ^' A . (5.10) 
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The corresponding variance 

/ + 00 
x^g^(x,t)dx^2Kj'', (5.11) 
-oo 

is characteristic of a general process of anomalous diffusion, precisely of slow diffusion for 
< a < 1, and fast diffusion for 1 < a < 2. 



The time-fractional diffusion equation. 

In literature there exist two forms of the time-fractional diffusion equation of a single 
order, one with Riemann-Liouvile derivative and one with Caputo derivative These forms 
are equivalent if we refer to the standard initial condition u{x,0) — uo{x), as shown in 
[47]. 

Taking a real number (3 e (0, 1), the time- fractional diffusion equation of order (3 in the 
Riemann-Liouville sense reads 



whereas in the Caputo sense reads 



^^tU = K,^, (5.13) 



where Kj^ is a sort of fractional diffusion coefficient of dimensions [Kp] — [L\'^[T]~^ — 
cm^/sec^. Like for diffusion equations of integer order (5.1) and (5.7)-(5.8), we consider 
the equivalent integral equation corresponding to our fractional diffusion equations (5.12)- 
(5.13), 

u{x, t) = uo{x) + j\t - rf-' ^^^^ dr . (5.14) 

The Green function ^^{x^t) for the equivalent Eqs. (5.12)-(5.14) can be expressed, also 
in this case, in terms of the M- Wright function, as shown in Appendix by adopting two 
different approaches, as follows: 



2 ^fWatf^l'^ ^' \ ^/K^t^/^ 



The corresponding variance can be promptly obtained from the general formula (5.5) for 
the absolute moment of the M- Wright function. In fact, using (5.5) and (5.15) and after 
an obvious change of variable, we obtain 

/+00 2 
^ x^ Gp{x, t) dx = YiJTi) ■ ^^-^^^ 
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As a consequence, for < ^5 < 1 the variance is consistent with a process of slow diffusion 
with similarity exponent H = /3/2. For further reading on time-fractional diffusion 
equations and their solutions the reader is referred, among others, to [39, 44, 45] and 
[62], [66]. 



The stretched time-fractional diffusion equation. 

In the fractional diffusion equation (5.12), let us stretch the time variable by replacing t 
with i"/^ where 0<Q;<2andO<^<l. We have 



namely 



(5-17) 



.D'-^— (5 18) 



where Ka/s is a sort of stretched diffusion coefficient of dimensions [K^p] = [L]^[T]~" = 
cm'^/sec°' that reduces to if /9 = 1 and to Kp ii a — /3. Integration of Eq. (5.18) gives 
the corresponding integral equation [57] 

nix, t) = Uo{x) + K^,^^^I^ ir/^ - r^/^-' dr , (5.19) 

whose Green function Qai3{x,t) is 



with variance 



/+00 2 
G^,p{x, t) dx = Y(j3Tr) ■ ^^-^^^ 

As a consequence, the resulting process turns out to be self-similar with Hurst exponent 
H — a/2 and a variance law consistent both with slow diffusion if < a < 1 and fast 
diffusion if 1 < a < 2. We note that the parameter ^ does explicitly enter in the variance 
law (5.21) only in the determination of the multiplicative constant. 

It is straightforward to note that the evolution equations of this process reduce to those 
for time- fractional diffusion if « = /3 < 1, for stretched diffusion if a 7^ 1 and /5 = 1, and 
finally to standard diffusion ii a = (5 = 1. 
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6 Fractional diffusion processes with stationary in- 
crements 

We have seen that any Green function associated to the diffusion-hke equations considered 
in the previous Section can be interpreted as the time-evolving one-point pdf of certain 
self-similar stochastic processes. However, in general, it is not possible to define a unique 
(self-similar) stochastic process because the determination of any multi-point probability 
distribution is required, see e.g. [58]. 

In other words, starting from a master equation which describes the dynamic evolution 
of a probability density function f{x,t), it is always possible to define an equivalence 
class of stochastic processes with the same marginal density function f{x,t). All these 
processes provide suitable stochastic representations for the starting equation. It is clear 
that additional requirements may be stated in order to uniquely select the probabilistic 
model. 

For instance, considering Eq. (5.18), the additional requirement of stationary increments, 
as shown by Mura et al., see [55, 56, 57, 58], can lead to a class {Ba,0{t), t > 0}, 
called "generalized" grey Brownian motion (ggBm), which, by construction, is made up 
of self-similar processes with stationary increments and Hurst exponent H = a/2. Thus 
{Ba,/3{t), t > 0} is a special class of H — sssi processes^, which provide non-Markovian 
stochastic models for anomalous diffusion, both of slow type (0 < a < 1) and fast type 
(1< a < 2). 

The ggBm includes some well known processes, so that it defines an interesting general 
theoretical framework. The fractional Brownian motion [fBm) appears for /3 = 1 and is 
associated with Eq. (5.7); the grey Brownian motion (gBm), defined by Schneider [64, 65], 
corresponds to the choice a — /3, with < ^ < 1, and is associated to Eqs. (5.12), (5.13) 
or (5.14); finally, the standard Brownian motion (Bm) is recovered by setting a = /3 = 1 
being associated to Eq. (5.1). We should note that only in the particular case of Bm the 
corresponding process is Markovian. 

In Figure 3 we present a diagram that allows to identify the elements of the ggBm 
class. The top region 1 < a < 2 corresponds to the domain of fast diffusion with long- 
range dependence^. In this domain the increments of the process Ba,p{t) are positively 

^According to a common terminology, H—sssi stands for iJ-self-similar-stationary-increments, see for 
details [70]. 

self-similar process with stationary increments is said to possess long-range dependence if the 
autocorrelation function of the increments tends to zero like a power function and such that it does not 
result integrable, see for details [70]. 
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correlated, so that the trajectories tend to be more regular {persistent) . It should be noted 
that long-range dependence is associated to a non-Markovian process which exhibits long- 
memory properties. The horizontal line a = 1 corresponds to processes with uncorrelated 
increments, which model various phenomena of normal diffusion. For a = (3 = 1 we 
recover the Gaussian process of the standard Brownian motion. The Gaussian process 
of the fractional Brownian motion is identified by the vertical line /3 = 1. The bottom 
region < a < 1 corresponds to the domain of slow diffusion. The increments of the 
corresponding process Ba^p(t) turn out to be negatively correlated and this implies that 
the trajectories are strongly irregular {anti-persistent motion); the increments form a 
stationary process which does not exhibit long-range dependence. Finally, the diagonal 
line {a = (3) represents the Schneider grey Brownian motion {gBm). 
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Figure 3: Parametric class of generalized grey Brownian motion 

Here we want to define the ggBm by making use of the Kolmogorov extension theorem 
and the properties of the M- Wright function. According to Mura and Pagnini [57], the 
generalized grey Brownian motion Ba^p{t) is a stochastic process defined in a certain 
probability space such that its finite-dimensional distributions are given by 

/,M:.,,:..,...,.„;7.,,) = ^===^y^ —,My.{—,)M,{r)ir, (6.1) 
with 

/ n \ 1/2 

i= \^T{1 + |3)-'Yl^^l^,f\t^^h)^3j . (6-2) 
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and covariance matrix 

laAt^tj) = + ^'^ ~ '^^ " ^^l"^' = 1' ■ ■ ■ (6-3) 

The covariance matrix (6.3) characterizes the typical dependence structure of a self-similar 
process with stationary increments and Hurst exponent H = a/2, see e.g. [70]. 

Using Eq. (4.18), for n = 1, Eq. (6.1) reduces to: 

1 1 

f^^^{x,t) = ^= / M^/2{\x\t-^/^r)M^{r,l)dr ^ -t--/^M^/2{\x\t-"/^) . (6.4) 
V At" Jo ^ 

This means that the marginal density function of the ggBm is indeed the fundamental 
solution (5.20) of Eqs. (5.17)-(5.18) with K^p = 1. Moreover, because Mi(r) = 6{t - 1), 
for P — 1, putting = 7a, we have that Eq. (6.1) provides the Gaussian distribution 
of the fractional Brownian motion, 

(2tt)-'^ If " V'^\ 

fa,l{Xi,X2, . . .,Xn;^a,l) = ^2 dct 7 ^^^^ ( ^ ^i^a^^i^h)^! j 1 ' (6-5) 

which finally reduces to the standard Gaussian distribution of Brownian motion as a = 1. 

By the definition used above, it is clear that, fixed /3, Ba^t) is characterized only by 
its covariance structure, as shown by Mura et al. [56], [57]. In other words, the ggBm, 
which is not Gaussian in general, is an example of a process defined only through its 
first and second moments, which indeed is a remarkable property of Gaussian processes. 
Consequently, the ggBm appears to be a direct generalization of Gaussian processes, in 
the same way as the M- Wright function is a generalization of the Gaussian function. 



7 Concluding discussion 

In this review paper we have surveyed a quite general approach to derive models for 
anomalous diffusion based on a family of time-fractional diffusion equations depending on 
two parameters a e (0, 2), /3 e (0, 1]. 

The unifying topic of this analysis is the so-called M- Wright function by which the 
fundamental solutions of these equations are expressed. Such function is shown to 
exhibit fundamental analytical properties that were properly used in recent papers for 
characterizing and simulating a general class of self-similar stochastic processes with 
stationary increments including fractional Brownian motion and grey Brownian motion. 
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In this respect, the M- Wright function emerges to be a natural generalization of the 
Gaussian density to model diffusion processes, covering both slow and fast anomalous 
diffusion and including non-Markovian property. In particular, it turns out to be the 
main function for the special H — sssi class of stochastic processes (which are self-similar 
with stationary increments) governed by a master equation of fractional type. 
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Appendix A: The fundamental solution of the time- 
fractional diffusion equation 

The fundamental solution ^^(x, t) for the time- fractional diffusion equation can be 
obtained by applying in sequence the Fourier and Laplace transforms to any form chosen 
among Eqs. (5.12)-(5.14) with the initial condition ^^(x,0+) = Uq{x) = S{x). Let 
us devote our attention to the integral form (5.14) using non-dimensional variables by 
setting = 1 and adopting the notation for the fractional integral. Then, our 
Cauchy problem reads 



g^{x,t) = 6{x) + J,^^{x,t). (Al) 

In the Fourier-Laplace domain, after applying formula (2.18) for the Laplace transform 
of the fractional integral and observing S{k) = 1, see e.g. [15], we get 

Gp{K,s) = ^Gis{k,s) , 

s s 



from which 



gp{K,s) = ^—^, 0</5<l, sft(s) >o, «:GIR. (A2) 

To determine the Green function Qii{x^t) in the space-time domain we can follow two 
alternative strategies related to the order in carrying out the inversions in (A. 2). 

(51) : invert the Fourier transform getting Qp{x, s) and then invert the remaining Laplace 
transform; 

(52) : invert the Laplace transform getting Gp{K, t) and then invert the remaining Fourier 
transform. 
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Strategy (SI): RecaUing the Fourier transform pair 



b + 26V2 
and setting a = 6 = s'^, we get 



^ ^e-l^l^'^', a,6>0, (A3) 



gp{x,s) = ]-s^l'-'e-\^\'^'\ (A4) 



2 

Strategy (S2): RecaUing the Laplace transform pair 

A Ep{-ct^), OO, (A5) 

sP -\- c 

and setting c = we have 

Gp{K,t) = Ep{-KH^). (A6) 

Both strategies lead to the result 



Gp{x,t) = ^M^/,{\x\,t) = \t-^'^Mp,^ 



\x\ 



(A.7) 



consistent with Eq. (5.15). Here we have used the M- Wright function, introduced in 
Section 4, and its properties related to the Laplace transform pair (4.15) for inverting 
(A. 4) and the Fourier transform pair (4.17) for inverting (A. 6). 



Appendix B: The fundamental solution of the time- 
fractional drift equation 

Let us finally note that the M- Wright function does appear also in the fundamental 
solution of the time-fractional drift equation. Writing this equation in non-dimensional 
form and adopting the Caputo derivative we have 

d 

*Dfu(x, t) — ——u(x, t) , — OO < X < +00 , t >0 , (B.l) 
ox 

where < /5 < 1 and u(,t,0^) = Uo{x). When Uo{x) = S{x) wc obtain the fundamental 
solution (Green function) that we denote by g^{x, t). Following the approach of Appendix 
A, we show that 

g;{x,t) = h'''^^{i^)^ iB.2) 

^ ( 0, X < 0, 

that for ^ = 1 reduces to the right running pulse 5{x — t) for x > 0. 
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In the Fourier-Laplace domain, after applying formula (2.19) for the Laplace transform 
of the Caputo fractional derivative and observing S{k,) = 1, see e.g. [15], we get 

/ G*^{k, s) - s^-^ = +mG*^{K, s) , 

from which 

g*U.s) = — , 0</3<l, > 0, K e IR . (B.3) 

' SP — IK 

Like in Appendix A, to determine the Green function t) in the space-time domain 

we can follow two alternative strategies related to the order in carrying out the inversions 
in (B.3). _ 

(51) : invert the Fourier transform getting Qjs^x, s) and then invert the remaining Laplace 
transform; 

(52) : invert the Laplace transform getting G*^{k, t) and then invert the remaining Fourier 
transform. 

Strategy (SI): RecaUing the Fourier transform pair 

^ ^e-^^ a,b>0,x>Q, (BA) 
h — in 

and setting a — s^~^, b — s^, we get 

^(x,.) = /-le-^^^ {B.5) 



Strategy (S2): Recalling the Laplace transform pair 

A Ep{-ct^), 00, (5.6) 

si^ -\- c 

and setting c = —ire, we have 

G}{K,t) ^ EpiiKt/^) . {B.7) 
Both strategies lead to the result (B.2). 

In view of Eq. (4.1) we also recall that the M- Wright function is related to the unilateral 
extremal stable density of index /3. Then, using our notation stated in [44] for stable 
densities, we write our Green function as 

g;{x,t) = ^x-'-'/n/{tx-'/^) , (5.8) 



To conclude this Appendix let us briefly discuss the above results in view of their relevance 
in fractional diffusion processes following the recent paper by Gorenflo and Mainardi [22] . 
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Equation (B.l) describes the evolving sojourn probability density of the positively oriented 
time-fractional drift process of a particle, starting in the origin at the instant zero. It has 
been derived in [22] as a properly scaled limit for the evolution of the counting number 
of the Mittag-LefHer renewal process (the fractional Poisson process). It can be given in 
several forms, and often it is cited as the subordinator (producing the operational time 
from the physical time) for space-time-fractional diffusion as in the form (B.8). For more 
details see [25], where simulations of space-time- fractional diffusion processes have been 
considered as composed by time-fractional and space-fractional diffusion processes. 

This analysis can be compared to that described with a different language in papers by 
Mccrschaert et al. [51, 52]. Recently, a more exhaustive analysis has been given by 
Gorenflo [17]. 
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